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ABSTRACT 

Dense stellar systems such as globular clusters, galactic nuclei and nuclear star clusters 
are ideal loci to study stellar dynamics due to the very high densities reached, usually 
a million times higher than in the solar neighborhood; they are unique laboratories 
to study processes related to relaxation. There are a number of different techniques 
to model the global evolution of such a system. We can roughly separate these ap- 
proaches into two major groups; the particle-based models, such as direct N— body 
and Monte Carlo models, and the statistical models, in which we describe a system 
of a very large number of stars through a one-particle phase space distribution func- 
tion. In this approach we assume that relaxation is the result of a large number of 
two-body gravitational encounters with a net local effect. We present two moment 
models that are based on the collisional Boltzmann equation. By taking moments of 
the Boltzmann equation one obtains an infinite set of differential moment equations 
where the equation for the moment of order n contains moments of order n + I. In our 
models we assume spherical symmetry but we do not require dynamical equilibrium. 
We truncate the infinite set of moment equations at order n — 4 for the first model 
and at order n — 5 for the second model. The collisional terms on the right-hand 
side of the moment equations account for two-body relaxation and are computed by 
means of the Rosenbluth potentials. We complete the set of moment equations with 
closure relations which constrain the degree of anisotropy of our model by expressing 
moments of order n + 1 by moments of order n. The accuracy of this approach re- 
lies on the number of moments included from the infinite series. Since both models 
include fourth order moments we can study mechanisms in more detail that increase 
or decrease the number of high velocity stars. The resulting model allows us to derive 
a velocity distribution function, with unprecedented accuracy, compared to previous 
moment models. 

Key words: stars: kinematics and dynamics - globular clusters: general - galaxies: 
nuclei 



1 MOTIVATION 

Statistical continuum models such as Fokker-Planck (FP) 
and moment models separate the treatment of the differ- 
ent astrophysical processes that control the evolution of 
the system. This allows us to isolate the effects of the dis- 
tinct dynamical mechanisms. In particular statistical mo- 
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ment models have provided us with important contribu- 
tions to the understanding of phenomena such as core col- 



lapse and gravothermal oscillations ( jBettwieser fc Sugimoto| 



19841. These models decompose the local velocity distribu- 



tion function into the different contributions of the moments, 
allowing us to study the influence of the different moments 
on the evolution of star clusters and the impact of different 
dynamical mechanisms on the moments of the distribution 
function. This has a bearing in a number of crucial problems 
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such as the contribution of high velocity stars to the evolu- 
tion of star clusters, which we only can address by including 
fourth order moments. 

We present in this paper two statistical moment mod- 
els for dense, non-rotational and spherically symmetric stel- 
lar systems, such as globular clusters (GCs) or nuclear star 
clusters (NCs). The models include fourth order moments 
and thus allow us to study astrophysical scenarios that af- 
fect the number of high velocity stars. The models describe 
the evolution of a stellar system that slowly evolves due 
to the effects of two-body relaxation. Moment models have 
the advantage over particle-based techniques in that they 
are computationally much cheaper, being based on the nu- 
merical integration of a relatively small set of partial dif- 
ferential equations with just one variable, the radius r. The 
numerical solution of the model equations is usually very 
fast as they are equivalent to one-dimensional hydrodynam- 
ical equations. Since the system is treated as a continuum, 
all macroscopic quantities (such as density, pressure and en- 
ergy flux) are smooth functions of radius r and time t and 
do not suffer from the characteristic noise of particle-based 
approaches. 

Moment models began with simple collisionless mod- 



els and progressed to the anisotropic gaseous model (Bet- 



twieser fc Spurzem 1986; |Louis fc Spur zem 199 1] |Spurzem 
1992||Giersz fc Spurzem|1994||Spurzem fc Takahashi|1995 l 



They have significantly contributed to the understanding of 
stellar dynamical systems by gradually adding new phenom- 
ena such as two-body relaxation, three-body encounters, and 
energy transport processes in stellar systems with a mass 
spectrum. 

Moment models could quite easily be coupled with hy- 
drodynamical solvers to simulate the dynamical evolution 



of dense gas-star systems (DGSS) in galactic nuclei (Lang- 


bein et al. ||1990| |Amaro-Scoane et al.|2001 


2002| |Amaro- 


Seoane & Spurzem|2004 Spurzem et al.|2004 


. In Langbein 


et al. |(|1990|) it was shown that gaseous models of dense star 



clusters can be regarded as a generalization of the Tolman- 
Oppenheimer Volkoff equation for relativistic anisotropic 
gases. Many years ago |Bisnovatyi-Kogan fc Sunyaev | ( |1972[ ), 



Vilkoviski ( 1975 1, and Hara ( 1978 1 have proposed DGSS as 



energy sources in galactic nuclei. Nowadays the idea is being 
reconsidered that supermassive stars are progenitors of the 
first supermassive black holes in galactic nuclei (Bcgclman 



20101, and that galactic nuclei in their variety of appear- 



ances could be determined by the interplay of stellar and 



gas dynamics, including star formation and feedback (Ciotti 
|et al |2009||Shin et al.|2010||Ciotti et al.|2010| . These topics 



deserve further investigation with improved stellar dynami- 
cal modelling, as we provide it here with our new momentum 
model. Therefore we think a fresh look at and improvement 
of the momentum model is timely and very useful. It should 
be noted that spherical symmetry yet has been a limitation 
of gaseous or momentum models of star clusters. However, 
also here a generalization at least to axisymmetric models is 
possible by describing viscosity through two-body relaxation 
in analogy to heat conduction (Goodman 1983a ( unpublished 
Ph.D. thesis). We have demonstrated that the aforemen- 
tioned Goodman models can be used and solved numerically 
with sufficient accuracy in the case of direct solutions of the 
orbit averaged Fokker-Planck equation ( jEinsel fc Spurzem| 
[19991 |Kim~et al. ||2002"1 |Kim et al.|[2004[ |Kim et al. ||2008| 



Fiestas fc Spurzem |2010| ). There is no reason to assume that 
also our momentum or gaseous model could not be extended 
to axial symmetry in the future, using appropriate implicit 
hydrodynamic solvers. 

By extending the model with additional equations cou- 
pled with collisional terms, we are in the position to address 
new problems. Thus we can investigate accretion theory 
(Amaro-Seoane et al. 2004[ ), stellar collision, gas dynam- 



ics and coupling with the stellar system, including radiative 



transfer and turbulences, the role of the loss-cone 


Amaro- 


Seoane et al. 2003 ; Amaro-Seoane & Spurzem 


2004 


Amaro- 


Seoane|2004 1 and tidal fields ( 


Spurzem et al. 


2005 


. Higher 



order moments are necessary to have a more realistic de- 
scription of the velocity distribution function and a more 
accurate description of relaxation, reducing the number of 
approximations necessary to the model. 

The numerical models used to study dynamical pro- 
cesses have to be constrained by comparison with observa- 
tions. In order to do so, both models and observations must 
fulfill certain accuracy requirements. There are many meth- 
ods for modeling GCs which can be separated into particle 
based methods such as ./V-body or Monte-Carlo simulations 
and continuum methods such as Fokker-Planck or moment 
models (see next section). In statistical moment models, we 
employ velocity moments to characterize the local velocity 
distribution function. The n-th moment of a velocity dis- 
tribution f(v) is defined as (v n ) = J v n f(v) dv (see also 



definition ( 14 1 ) . The accuracy of these models is then lim- 
ited by the order of the highest moment included to describe 
the velocity distribution. A physical interpretation for each 
moment up to the fourth order can be given. Since each 
stellar dynamical process driving the evolution of a cluster 
has a different impact on the local velocity distribution, this 
motivates us to construct a distribution function that is able 
to reflect the effects of each of these processes properly so 
as not to lose information that influences the clusters evo- 
lution. The velocity distribution can be written as a series 



expansion using a truncated Gauss-Hermite series ( Gerhard 
1993 van der Marel fc Franx|1993 1 to illustrate the meaning 
of the first four moments: 



f(ur) oc exp(- 



2(7 



(1) 



v r might be the velocity in radial direction (or the line-of- 
sight velocity which is the velocity measured in direction of 
an observer). v r , a, /13 and /14 are free parameters and will 
be explained in the following. 

• 0th moment: 

The zeroth moment of a velocity distribution is 1 due to 
normalization. 

• 1st moment: 

The first moment of a velocity distribution is the mean ve- 
locity v r and denotes the bulk mass transport velocity. 

• 2nd moment: 

The second moment of a velocity distribution is the variance 
o and is equal to the velocity dispersion. It determines the 
width of f(v T ) and thus the scattering of stellar velocities 
around the mean velocity v r . If f(v r ) is fully determined 
by v r and a and /13 = /14 = it is a Gaussian (top pannel 
in figure [TJ corresponding to thermal equilibrium. Then the 



Higher order moment models of dense stellar systems 3 



symmetry of the one-dimensional velocity distribution f(v r ) 
to v r reflects isotropy. 

• 3rd moment: 

The third moment, denotes the transport of random kinetic 
energy and depends on /13 . If the third moment of the veloc- 
ity distribution does not vanish, implying that /13 7^ 0, then 
the shape of the velocity distribution is a skewed Gaussian 
(figure [l] upper middle pannel) . The asymmetry indicates 
the direction of the energy flux, and the uneven distribution 
of velocities in different directions denotes anisotropy. 

• J^th moment: 

The fourth moment is a measure of the excess or deficiency 
of particles/stars with high velocities as compared to ther- 
modynamical equilibrium, and depends on the value of /14. 
An excess of particles with high velocities results in thicker 
wings of the velocity distribution and a more pointed maxi- 
mum (figure[lj lower middle pannel). A deficiency of high ve- 
locities causes a broader shape around the mean and thinner 
wings of the velocity distribution (figure [T] bottom pannel). 

Third and fourth order moments therefore denote de- 
viations from thermodynamical equilibrium. Modeling pro- 
cesses that lead to the transport of random kinetic energy in 
a cluster or that strongly affect the high velocity wings of the 
distribution suggest the use of a model that includes 4th or- 
der moments. These processes are, for example, the "evapo- 
ration" of high velocity stars from the cluster, which reduces 
the number of high velocity stars. On the other hand, bi- 
naries and a mass spectrum transfer kinetic energy between 
different stellar components and thereby produce high veloc- 
ity stars. These high velocity stars then transfer their excess 
energy to their environment in subsequent distant two-body 
encounters which can lead to a transport of kinetic energy 
between different regions in the GC. 

Neglecting third and fourth order moments in these 
cases results in a loss of information by failing to fully model 
the effect of the processes they represent on the evolution of 
the cluster. 



2 PARTICLE-BASED TECHNIQUES VS 
STATISTICAL METHODS 

The methods for studying star clusters can be divided into 
two types; statistical continuum models, such as Fokker- 
Planck (FP), or moment models and particle-based tech- 
niques, such as direct iV-body models and Monte Carlo. 
They have different advantages and deliver complementary 
information about the processes and mechanisms that drive 
the evolution of star clusters. 

2.1 Direct integration techniques 

By using direct iV-body we integrate Newton's equations of 
motion. In principle all gravitational dynamics phenomena 
are naturally included in the integration. Thus, this method 
is not subject to any approximations nor restricted to any as- 
sumptions, such as spherical symmetry. In contrast to statis- 
tical methods, it does not require additional physics in order 
to include gravitational interactions between pairs, triples 
(binary-star interactions) or quadruples (binary-binary in- 
teractions) as they are inherent to the model. Including a 
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Figure 1. These four plots show one-dimensional velocity distri- 
bution functions for different cases, top: Gaussian velocity distri- 
bution describing thermodynamical equilibrium with a variance 
of <t — 10km/s. The Gaussian appears in the subsequent panels 
for comparison (black), upper middle: velocity distribution (grey) 
with a skewness in positive u r -direction indicating energy flow in 
rv-direction. lower middle and bottom: two velocity distributions 
(grey) with an excess and deficit of high velocity stars respectively 
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mass spectrum or tidal field is also, in principle, straightfor- 
ward. On the other hand, direct-summation methods of this 
type are computationally expensive, and as a consequence it 
is not possible to realistically model a stellar cluster with a 
typical number of 10 7-8 stars. This is due to the fact that the 
computation of all pairwise interactions of a system consist- 
ing of TV particles scales with TV 2-3 . Using modern hardware 
we are severely limited to integrations of at most a few 10 6 
particles for a very short time, typically a few dynamical 
times. Another drawback of direct TV-body is that it suffers 
from noise, as an individual N-body calculation in star clus- 
ter dynamics have exponential instabilities; nevertheless, the 
results can be used in a statistical average (e.g. |Miller|T964[ 
Giersz fc Heggie|1994a l. 

There exist many schemes for integrating Newton's 
gravitational equation, some of them are faster and more 
effective than others. Among these we should mention the 
Euler scheme or an improvement of this, the leapfrog scheme 
(e.g. Hut et al.|1995 1. We can gain more accuracy by the 



( 1971a|b 1972 19751 and later improved by Stodolkiewicz 



vided difference scheme or the Hermite scheme (Makino & 
Aarseth|199"2] |Aarseth|T999] >, which is used in the NBODY6 



and NBODY6++ codes for the orbit integration. Addition- 
ally, various TV-body codes incorporate a number of ap- 
proaches which are necessary for maintaining adequate ac- 
curacy and efficiency over many dynamical times; these in- 
clude the use of many individual time steps, computation 
of forces from near neighbors and distant stars with differ- 
ent frequencies, special treatments of compact pairs (bina- 



ries) and other few-body configurations ( Mikkola & Aarseth 



1990 1993 1. Direct TV-body simulation is a powerful tool for 



realistically simulating a wide range of astrophysically in- 
teresting scenarios such as black holes in galactic nuclei or 
GCs, binaries of massive black holes in (rotating) clusters 
(Amaro-Seoane & Freitag 2006 Amaro-Seoane et al.||2009 
Amaro-Seoane et al.||2010| or binary black hole mergers in 



galactic nuclei ( jBerentzen et al.||2009[ ) 



(1982 [19861. In contrast to the models of Spitzer, Henon's 
models assumed dynamical equilibrium; the distribution 
function must also depend only on isolated integrals of mo- 
tion. It is worth mentioning that it was the first scheme 



to break through the impasse of core collapse ( Henon 



The algorithm was further improved by Stodolkiewicz 



1975) 



( 1985 1 



by including processes such as the formation of binaries by 
two- and three-body encounters, mass loss from stellar evo- 
lution and tidal shocking. 



Giersz ( 1998 



2001 \ in a series of papers modeled uCen 
( Giersz fe Heggjep003| , M4 ( |Heggie fc Glers^l[2008| M67 
( Giersz et al.|2008[ ) and NGC 6397 ( |Giersz fc Heggiel[2009| 



with MC techniques. In these papers several additional im- 
provements were also included, as two-body relaxation, most 
kinds of three- and four-body interactions involving of pri- 
mordial binaries and those formed dynamically, the Galac- 
tic tide and the internal evolution of both single and binary 
stars. MC techniques can be coupled with continuum models 
to describe the stochastic process of binary formation energy 
generation and movement (Spurzem & Giersz||1996 Giersz 



& Spurzem||2000 2003 1 . This has been successfully used to 



examine the gravitational radiation from binary black holes 
in star clusters ( Downing et al.||2009 l 



Joshi et al. (2000 2001 1; Fregeau et al. (20031; Fregeau 



& Rasio ( 2007 1 developed a MC technique based on a mod- 



ified version of Henon's algorithm for solving the Fokker- 
Planck equation. Their scheme includes a mass spectrum, 
stellar evolution, and primordial binary interactions and 
the direct integration of binary scattering interactions. The 
Henon-type MC approach has been used by M. Freitag, who 
developed another MC code with the special purpose of 
studying semi-Keplerian systems. Applying this code he ex- 
tensively studied the structure of galactic nuclei containing 
a central M BH ( |Freitag|200"o"] |Freitag fc Benz|2002| |Freitag 
|et al.|2006| |. 



2.2 The Monte Carlo approach 

Other powerful particle-based techniques are the Monte 
Carlo (MC) methods, in which relaxation is treated using 
the Fokker-Planck approximation. These methods rely also 
on the assumptions that the system is spherically symmetric 
and that the gravitational potential can be separated into 
two parts. The advantage of MC is that it is orders of mag- 
nitude faster than direct TV— body, yet it is still slower than 
statistical methods and also suffers from numerical noise. 
Spitzer and collaborators pioneered the MC scheme in a 



series of papers, such as|Spitzer & Hart|i 


1971a|b|);|Spitzer& 


Shapiro (19721; Spitzer & Thuan (1972 


; Spitzer & Cheva- 


lier (19731; Spitzer & Shull (1975a|bl; 


Spitzer & Mathieu 



and his collaborat ors ([Shapiro fc Marchant|1978||Marchant| 
fc Shapiro| [19791 |1980| |Duncan fc Shapiro||1982| |Shapiro| 
1984[ ). MC, being particle-based, follows the individual stel- 



lar orbits and allows us to model processes occurring on 
both relaxation and crossing time scales. Spitzer's method 
was used to explore a variety of important phenomena, in- 
cluding mass segregation, anisotropy of the velocity distribu- 
tion, tidal shocking, and the role of primordial binary stars, 
to mention a few. 

The second MC approach was devised by |Henon| 



2.3 A statistical model: 
technique 



The Fokker-Planck 



Fokker-Planck (FP) models are based on the direct numeri- 



cal solution of the orbit-averaged FP equation. Cohn ( 1979 
1980) pioneered a direct numerical finite-difference solution 
of the 1-dimensional FP equation (for a phase space distri- 
bution function: / = f(E)). Similar methods had been de- 
eloped for a fixed potential by |Ipser| ( |1977[ ) and by |Cohn"fc 



Kulsrud ( 1978 1, and since then different FP codes have been 



written independently by Inagaki & Wiyanto ( 1984 1 and by 
|Chernoff fc Weinberg! ( |1990[ ). Whereas Cohn's formulation 
assumes spherical symmetry, codes which can handle a ro- 
tating cluster have been devised by Goodman ( 1983a I and 
|Einsel fc Spurzem| ( p96] . |Takahashi| ( |1995[|1996||1997| ) has 



developed FP models for GCs, based on the numerical solu- 
tion of the orbit averaged 2D FP equation (i.e. solving the 
FP equation for the distribution / = f(E, J 2 )) as a function 
of energy and angular momentum, and thus accounting for 
anisotropy. 

|Drukier et al.| ( |1999[ ) followed with results from another 
2D FP code based on the original idea of Cohn ( 1979 1 
There have been several comparati ve studies (|Giersz fc Heg 



gie|1994a|b|[Wn|Taka"hashi|1995||Giersz fc Spurzem|1994 



Spurzem fc Takahashi|1995[|Freitag et al.|2006||Khalisi et aL] 
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20071 showing that for isolated, non-rotating star clusters 
the results of FP simulations are generally in good agree- 
ment with those of iV-body simulations. However, when a 
tidal boundary is included, discrepancies between iV-body 
and FP models occur. 



Also, Einsel & Spurzem ( 1999) found that rotating GCs 



collapse faster than non-rotating ones with a 2D FP tech- 
nique that had a distribution function depending on the z- 
component of the angular momentum, / = f(E,J z ). Kim 
( 2002 ) improved the approach by including an energy 



et al. 



source due to formation and hardening of three-body bina- 
ries. These two studies only investigated single-mass models. 
Later, Kim et al. ( 20041) extended this method to multi-mass 



systems, finding interesting results concerning segregation of 
mass and angular velocity with heavy stars in the cluster. 
Fiestas et al. ( 2006 1 have modeled rotating globular clus- 



ters and Fiestas & Spurzem (2010) included a star accreting 



black hole with a loss cone. Comparative studies for rotat- 
ing star clusters between FP and A-body methods have been 



done as well by Boily (20001; Boily & Spurzem 



et al. (20051; Ernst et al. (20071; Kim et al. (20081. They 



(2000); Ardi 



produced fairly similar results, though there were small dis- 
crepancies in the core-collapse time. 

2.4 Advantages and disadvantages of statistical 
models as compared to direct-summation 
techniques 

/.From the three different techniques, direct iV-body models 
appear as the most realistic model. However, as mentioned 
before, it suffers from exponential instabilities; small devi- 
ations in the initial conditions result in exponential diver- 
gence of the phase space distribution of the particles of the 
system ( |Miller|[T964] |Giersz fc Heggie||1994a| ). These insta- 
bilities make it difficult to compare a realistic model of GCs 
to observational data. Statistical models produce averaged 
physical quantities and are better suited for comparison with 
observations. 

Also, as A-body models are not restricted by boundary 
conditions such as spherical symmetry, they can be applied 
to the widest range of stellar dynamical systems and study 
them under the most diverse scenarios. This relies on a mi- 
croscopic description of dynamical processes and translates 
into a complexity that requires a massive computational ef- 
fort. As a consequence we depend on the development of 
hardware to push the number of particles that we can in- 
tegrate forward. On the other hand, statistical continuum 
models which are based on a comparatively small set of dif- 
ferential equations are computationally cheap. 

These algorithms have also the important property that 
the contribution of various dynamical processes to the over- 
all evolution of a star cluster can be isolated. This is so 
because the different mechanisms have to be included sep- 
arately by additional terms in the model equations. There- 
fore, it is possible to identify each mechanism and its effect. 

The downside of statistical moment models is that they 
are subject to a large number of approximations. Some of 
these approximations are inherent in the approach, such as 
the description of the phase space distribution function by 
a finite number of its velocity moments. Additional ap- 
proximations consist of the limitation in the number of pro- 
cesses included such as two-body relaxation, star-binary de- 



flections, binary-binary encounters or anisotropy. Such pro- 
cesses are natural in A-body models. The bottom line is 
that in order to build up a detailed understanding of stel- 
lar dynamical systems we need the different properties of 
particle-based and statistical models. 



3 SELF-GRAVITATING, CONDUCTING GAS 
SPHERES 

In the previous section we have given an overview on the 
different numerical tools to address stellar dynamics includ- 
ing relaxation. Now that we have highlighted the advan- 
tages of statistical methods, we introduce an interesting al- 
ternative to FP . More than 35 years ago |Hachisu et al.| 
( |1978[ ) and |Lynden-Bell fc Eggleton| ( [l980| proposed trans- 



port process in a self-gravitating, conducting gas sphere as 



twieser (19831; Bettwieser & Sugimoto (19841; Bettwieser 


& Spurzem (19861; Heggie ( 


19841 Heggie & N. Ramamani 


(|1989|) and|Louis & Spurzem 


(|1991|) implemented anisotropy 



and Giersz & Spurzem ( 1994 1 and Spurzem & Takahashi 
( 1995[ ) added a multi-mass distribution and improved the de- 



tailed form of the conductivities to have better accuracy. The 
resulting model is often called the anisotropic gaseous model 
(AGM). This allows us to compare with N— body models 
to calibrate the approach. |Amaro-Seoane et al.| ( |2004[ ) ad- 
dressed the accretion of stars on to a massive black hole by 
adding collisional terms corresponding to loss-cone physics 
as well as tidal effects and |Spurzem et al.| ( |2005| ) investi- 
gated the evolution and dissolution of star clusters under 
the combined influence of internal relaxation and external 
tidal fields. 

In this approach, we emulate spherically symmetric sys- 
tems as a continuum; relaxation is treated as a diffusive pro- 
cess in phase space using the FP equation. We employ the 
local approximation to simplify the FP equation by neglect- 
ing the diffusion in position. The idea behind this is that 
an encounter takes place in a volume that is much smaller 
than the dimensions of the whole system. We model energy 
transfer by a local heat flux equation with an appropriately 
tailored conductivity. 

The basis of the equations of the model is the FP equa- 
tion which describes the time evolution of the probabil- 
ity density function. Using spherical polar coordinates, the 
Boltzmann equation takes the form: 



9f , df df df df 

-TTT +V T —- +V T h V$ h V S 

at or ov T 0V4, ove 



(2) 



In the last equation the right-hand side denotes that colli- 
sions are given in terms of the FP approximation. Due to 
symmetry, we can define a tangential velocity v\ — vl + Vg , 
so that we have two velocities v t and v r to describe the sys- 
tem. The "centralized" moments are defined by multiplying 
the velocity distribution / with powers of v t and (v r — v r ) 
and integrating over velocity space. 

The term "centralized" means that the moments are 
defined with respect to the mean velocity components v r — 
(v r ) = u and vt = (v t ) = 0, because we assume spherical 
symmetry. The order of a moment is defined by n+m where 
n and m are the powers of velocities in the definition of mo- 
ments, i.e. J (v r — v r ) n v^ f d 3 v. The moments defined this 
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way correspond to the density of stars, p, the bulk velocity, 
u, the radial and tangential pressures, p r and p t , and the 
radial and tangential kinetic energy fluxes, F r and F t . In 
order to obtain the set of differential moment equations, we 
multiply equation (J2J with powers of Vt and (v r — v r ) and 
integrate it in velocity space. After some recasting the inte- 
grals can be substituted by the moments. Up to second order 
the moment equations are the continuity equation, the Euler 
equation (force) and radial and tangential energy equations: 



dp 

at 
du 
at 

dp r 
dt 



+ 



1 d_ 

dr 



(r up) = 



du GM r 1 dp r p r — Pt 



+ 



dr 
1 d 



+ 



+ 



p dr 



+ 2'- 



.2 Q r 

-Pt 



(r up r 



+ 2p, 



du 
dr 



+ 



pr 
1 d_ 
dr 



(r 2 F r 



2F t 



+ 



.2 Q r 
Pi 



+ 



Spr 
St 



(3) 



+ 



(r up t ) + 

Spt 
It 



bin3 

2p r u 
r 



2r 2 dr (V *■ 



_3p r 

5" 

dpt 

~&t 

3 Pr 
10 X A tr 

Here Xa is a numerical constant related to the time-scale 
of collisional anisotropy decay, necessary to describe the re- 
laxation effects on cluster evolution. It should become unity 
when describing GCs using higher moment models. How- 
ever, this can only be confirmed by simulations. The value 
of Xa is discussed in Giersz Spurzem] ( |1994[ ) and is cho- 
sen by calibrating with N— body simulations. The authors 
found that Xa = 0.1 is a physically realistic value within the 
half-mass radius for all numbers of particles. 

The two terms on the right-hand sides of the equations 
for radial and tangential energy equations are the collisional 
terms. The fist term accounts for relaxation from uncorre- 
cted two-body encounters and can be derived from the FP 
equation. The second term, which is marked with "bin3", 
refers to star-binary encounters. Close 3-body or star-binary 
encounters generate kinetic energy. If the energy generation 
is high enough, this mechanism can reverse core collapse. 

The radial and tangential pressure, p r and p t are re- 
lated to the random velocity dispersions; p r = po 2 and 



more equations are set up, a mass relation and two equations 
accounting for heat flux. The mass relation defines the mass 
M r contained in a sphere of radius r, 



dM r . 2 

— — = 47IT pM 

dr 



(5) 



where pM — M ■ p is the mass density and M the mass of 
the stellar component. We thus obtain a set of gas dynam- 
ical equations |3| coupled with the Poisson equation (25 1. 
Since the moment equations of order n obtained from the 
Boltzmann equation contain moments of the order n + 1, 
we need closure relations connecting the moments of order 
n + 1 with lower-order moments. This is achieved with the 
heat conduction closure, a phenomenological approach ob- 
tained in an analogous way to gas dynamics. It is motivated 
by the resemblance between a star consisting of a large num- 
ber of atoms and a star cluster with large number of stars 
not only on the simple level of the virial theorem but also 
due to similarities in heat transport, energy generation and 
core-halo evolution. It was used by Lyndcn-Bcll & Eggleton 
(19801, initially restricted to isotropic systems. In this ap- 



proximation we assume that heat transport is proportional 
to the temperature gradient 



F = - k jt = - A ~<r- 

dr dr 



(6) 



This equation describes the heat flux in gases and liq- 
uids and for this reason the models using this closure are 
also called conducting gas sphere models. Even though the 
use of equation (J6j) is based on the assumption of small mean 
free paths for the particles, which is certainly questionable 
for stellar dynamical systems, models like the AGM agree 
with other modeling methods (e.g. iV-Body, FP) ( | Giersz fc| 
Spurzeni|T9 94 Spur zem fc Takahashi|1995[ ) 

In the classical approach A oc pX /t, where A is the 
mean free path and r the collisional time. Choosing the 
Jeans length X 2 j — a 2 / (AivGp) for A 2 and the standard Chan- 
drasekhar local relaxation time t rx oc <r 3 /p ( Chandrasekhar 



1942[ ) for r, we obtain the conductivity A oc p/a. More pre- 
cisely, the conductivity takes the form found in |Lynden-Bell| 
fc Eggleton|(|1980||: 



lar clusters such as the radial velocity dispersion. The av- 
erage velocity dispersion is a 2 = (a 2 + 2cr 2 )/3, where the 
factor 2 comes from the fact that there are two tangential 
directions. The radial energy flux of random kinetic energy 
is F — (F r + Ft) 12. We can see this by adding the two- 
moment equations for radial and tangential pressure to ob- 
tain the gas-dynamical equation for the energy density. The 
velocities for energy transport are defined by 



v t = 



Fr_ 
3p r 

2p, 



+ u 



+ u 



(4) 



In the case of we 



isotropy, p r — pt everywhere and 
hence F r = t-Ft, so that v r — v t . Therefore, the transport 
velocities for radial and tangential random kinetic energy 
are equal. 

In order to close the set of moment equations (|3| three 



A = 



3CGmpN 



(7) 



where C is a dimensionless numerical constant of order unity. 
By means of the velocities of energy transport the heat flux 
equation can be recast to find the two closure relations in 
the anisotropic case 



U + 



X da 2 
4nGpt rx dr 



= 



Vr =Vt, 



where 



A : 



27V7T 

10 



c 



(8) 



(9) 



It should be emphasized that A is a free parameter 
that has to be determined by comparison with other mod- 
els such as iV-bod y (jGiersz fc Spu rzem 1994]), Louis' fluid 
dynamical model (Louis & Spurzem 19911 or FP models. 
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In the isotropic limit, A is just a scaling scaling factor, but 
when taking into account anisotropy, A prescribes the rela- 
tive speed of two processes: the decay of anisotropy and the 
heat flow between warm and cold regions. With increasing 
A heat flows faster, so there is less time for gravitational 
encounters to destroy anisotropy. A larger A thus results in 
stronger anisotropy. 



4 HIGHER MOMENT MODELS 

In this section we present a new higher order moment model. 
We derive the model equations which consist of differential 
equations for the velocity moments of the phase space dis- 
tribution function, a Poisson equation and three equations 
to close the system of equations. We first compute the left- 
hand sides of the differential moment equation and then use 
a polynomial ansatz for the phase space distribution func- 
tion to obtain the right-hand sides. We define two models, 
model a and model b, which differ in the number of differ- 
ential (moment) equations and their closure relation. 



4.1 Left-hand sides 

Without collisions, the Boltzmann equation takes the form 
of a conservation equation (df /dt = 0) and describes the ad- 
vective rate of change of the phase space distribution func- 
tion /. If we follow the trajectory of a particle in a system de- 
scribed by the collisionless Boltzmann equation, the number 
density in phase space around the particle does not change. 
This implies that flow in phase space is incompressible. It 
becomes compressible when collisions are introduced with 
FP terms on the right hand side of the Boltzmann equation. 

Assuming that the stellar system is spherically sym- 
metric, we can use spherical coordinates when we write the 
collisional Boltzmann equation, 



df df . df . df . df 

dt dr dv r dvg ov$ 



(10) 



Using the Lagrangian of a particle in a spherical sym- 
metric potential $(?-, t), we have that 



+ r 2 sin 0</> ) - $(r,t) 



(11) 



We then apply the Euler-Lagrange equations to the La- 
grangian, to derive the equations of motion 



d$ Vg + v% 
dr r 



V<j, 



V r Vg 

r 



+ 



u 4> 
rta.n9 

vgv,/, 

rtaxiO 



(12) 



After substituting equation ( 12 1 into the Boltzmann 
equation ( | 10| , we use the approach of spherical symmetry 

to define the tangential velocity v t 



v f> + v % anc l obtain 



dt dr 



4 



9<f>\ df v r v t df 
dr J dv r r dvt 



(13) 



We now define the velocity moments of the distribution 
function / = f(r,v r , Vt,t) by multiplying it by powers of v r 
and vt and integrating over velocity space, 



n,m] = / d «/d™ v™ = 2ir 



rOG re 
JO J -t 



j p n 7 



n. m+1 



(14) 

Again, the order of a moment is defined as k = n + m. 
To obtain the differential equations for the moments [n, m] 
we multiply equation ( |13| l with powers of v r and v t and 
integrate over velocity space. After some recasting we can 
substitute the integrals by [n, m] which yields 



3, . d , , . m + 2 r . . 

— [n,m\ + — [n + l,m J H [n + l,m\ 

dt dr r 



n d<& 
~\n — 1, m + 2] + n\n — 1, ml — — 
r dr 



5 . 

-[n, m 



(15) 

We now want to find a differential equation equivalent 



to equation ( 15 \ for centralized moments. The centralized 
velocity moments are defined with respect to their mean 
velocity. Due to the assumed spherical symmetry of the 
system the mean velocities of the tangential components 
vg = v<f, — vanish. The mean velocity is only given by 
the radial velocity component 

Vr — [1)0] — 11 — 27T / dv t / dv r fv r Vt (16) 
JO J -oo 

We hence obtain the definition for centralized moments by 



substituting v r in equation (141 with (v r — v r ). Furthermore, 



the centralized moments can be expressed in terms of the 
moments [n, m] and are defined as 



(n, m) 



d v (v r ~ v r ) n v™ f 

dv r (V r — Vr)™ V^ +1 f (17) 

= £(?) (-I)*-* [l,0] n - fc [fc,m] 



2tt / dv t 



It is evident from the second line of (171 that the first cen- 
tralized moment (1,0) = 0. 

We adopt the following notation for the centralized mo- 
ments: 



p 


= (0,0) 


(1,0) 


= 




Vr 


= (2,0) 


2p t = 


(0,2) 




F r 


= (3,0) 


F t = 


(1,2) 






= (4,0) 


K r t = 


(2,2) 


K t = (0,4) 


Ct r 


= (5,0) 


G r t = 


(3,2) 


Gt = (M) 


H r 


= (6,0) 


H r ,t - 


= (4,2) 


Ht.r = (2, 4) 



(18) 

Again, p is the particle density, p r and pt are the radial 
and tangential pressure and are related to the radial and 
tangential velocity dispersion a r — p r /p and at = Pt/p, and 
F r and Ft denote the radial and tangential energy flux. 

We obtain a linear system of equations which can be 
solved for the moments [n, m] by computing all centralized 



moments (n, m) up to order n + m = 6 using equation ( 17 1 
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[2,0] 
[0,2] 
[3,0] 
[1,2] 
[4,0] 
[2,2] 
[0,4] 
[5,0] 
[3,2] 
[1,4] 
[6,0] 
[4,2] 
[2,4] 
[0,6] 



= p r + pU 

= 2pt 

= pU 3 + 3lip r + F r 

= 2up t + Ft 

= pit 4 + 6u 2 p r + 4uF r + K T 
= 2u 2 p t + 2uF t + K r t 
= Kt 

= pU 5 + Wu S p r + 10lt 2 -F r + 5uK r + G r 

= 2u A p t + 3u 2 F t + 3uK r t + Grt 

— UKt + Gt 

— pu 6 + 15u 4 p r + 20u 3 F r + 15u 2 ft r + QuG r 

= 2u i p t + ^U 3 F t + §UK T t + 4uG r t + H rt 

— U K t + 2uGt + H tr 

= H t 



(19) 



To obtain the differential equations for the centralized 
moments, we substitute the transformation from equation 



(191 into equation (151 and then successively use differential 



equations for lower-order moments to simplify the differen- 
tial equations for higher orders. We divide the differential 
moment equations into three sets, defined as follows 



Set I 
dp 



dt 



+ div(ptt) = 



5_p 

St 

dp 



^ + div(p V) + — ^ 
+ div(F r + up r ) + 2p, 



du 
dr 



9$ 



Spu 



-F t 



dt 

2% + divfFi + 2up t ) + -(F t + 2up t 
at r 



Sp r 
St 



Spt 
St 



(20) 



Set II 

OF, 



du 



, + div(«v + uF r ) + 3F r „ 
at or 



3— divp r 



2p r p t , _ / SFr 

~^>-y~5t 



9F t .. , _ . _ du 2p t .. 

— - + dav(Krt + uF t ) + F t - divp r 

at or p 



~dt 



- (K t — 2Krt — 2uF t 



+ div(G r + 



4 (G r t 



2p t F r 



■ 4k, 

IT 



du 
dr 



8 Ft 
St 

F r 

4 — divp r 
P 



(21) 



+ div(G r t + uK rt ) + 2n rt ^ 2 — divp r 

at dr p 



+ -(G rt -Gt+UKrt+2 — ) 

r p 

^ +div(G t +UK t ) + 4 (G t +uK t ) 
dt r 



8 Krt 

~ST 



enc 
8Kt 
lit 



Set III 

dG 



)t +div(H r +uGr)+5Gr^r- 

dt dr 



5 — divp r 
P 



r 

dGrt 



(Hrt 



PtK r , _ ( 5G r 



(22) 



+ div( H rt + uGrt) + 3G rt ^r ~ 3— divp,. 
dt dr p 

+ ~(2Hrt - 3Htr + 2uGrt + & P —) = ( ^ 

r p \ dt 

dGt , TT „ , „ du n t ,. 

— h div(_H tr + uGt) + Gt— divp r 

dt dr p 

_I (jfft _ 4jfftr _ 4u G t -2^)=(^ 
r p \ St 

Note that the divergence operator in spherical symme- 
try reduces to: 



,. 1 d 2 

div = — — r 
r 4 dr 



(23) 



We now define two models with different accuracies, 



Model a - including ( 20 1 and 



Model b - including all, (201, (211 and (22 I 



21) 



The potential $(r, t) is determined by the fraction of 
cluster mass M r (t) contained at radius r 



$ = 



GM r 

r 



(24) 



$ obeys the Poisson equation A<E> = 4npM, where pm = Mp 
is the mass density of the cluster. This leads to the equation 
for M r 



dM r . 2 

— — = 4nr p M 
dr 



(25) 



We note that the moment equations of order n contain 
moments of order n + 1. To close the system of equations 
we need closure equation where moments of order n + 1 
are expressed with lower-order moments. We derive these 
relations in the next section. 



4.2 FP collision terms 

We now compute the right-hand sides of the differential mo- 
ment equations ( |20[ ), | |21[ ) and ( |22[ ), i.e. the collisional terms. 
Our starting point is the collisional Boltzmann equation 
1 10 1. We have to find an expression for the term (8f/St) cnc . 
This can be done by approximating it with the Fokker- 
Planck equation, which requires that the evolution of the 
stellar system is driven by uncorrelated distant encounters. 

The Fokker-Planck equation is a diffusion equation that 
describes the diffusion of the phase space distribution func- 
tion in position and velocity space. We assume that the vol- 
ume in which a stellar encounter takes place is small when 
compared to the volume of the whole system. As a conse- 
quence, we can assume that during an encounter only the 
velocity of the particle is modified, but not the position. We 
thus neglect the diffusion of the phase space distribution 
function in position space. This approach is usually referred 
to as the "local approximation". Therefore, the right-hand 



side of equation (101 is 
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'ST 
St , 



= -E 



8V; 



(f(x,v)D(A Vi )) 



+ 9 E 



dvidi 



(f(%, v)D(AviAvj)) 



(26) 



D(Avi) and D(AviAvj) are the diflusion coefficients 
which depend on position and velocity coordinates. They 
determine the diffusion of the phase space distribution func- 
tion in velocity space and describe the average change of 
the i-th component of velocity per unit time due to stel- 
lar collisions. This is expressed by their dependence on the 
change of the ith velocity component Avi . Note that there 
are no diffusion coefficients that depend on Axi as we are 
using the local approximation. The diffusion coefficients are 



(Rosenbluth et al. 19571 



D(Avi) = 47vG 2 m f \nA^-h(v) 



D(AviAvi) = 4nG mf In A 



dvi 
d 2 



(27) 



dvidvj 

In A is the Coulomb logarithm, where A is the ratio between 
the upper and lower cut-off impact parameter 6 in a stellar 
collision. h(y) and g(v) are called the Rosenbluth potentials 
which are given by 



h(v) = (m + to/) 



\v - Vf\ 



g{y) = mj j f(v~})\v~Vf\d 3 v f 



(28) 



Thus, to denotes the mass of a test star that moves through 
a distribution f(Vf ,(!,/) of field stars with a mass to/. The 
FP equation then takes the form: 



- v r — V cos 6' 
ve — V sin 6 cos cj>' 

V,j, 



(30) 



V sin 6' sine//, 



where V is the modulus of v , 0' the angle between v and the 
radial direction, and <j>' the angle which defines the orienta- 
tion of the tangential component of v in the vg-y^-planc. 
Note that our model describes non-rotational spherically 
symmetric systems and thus the mean tangential velocities 
are vg — v$ — 0. Thus, equations (301 denote the compo- 



nents of the radial and tangential velocities with respect to 
their means. As the phase space distribution function / only 
depends on Vg + iA, it is independent on the angle <j>' . We 
can henceforth omit the prime from angle 9' . 

Since we are operating in velocity space, in the following 
we refer to / as the (local) velocity distribution function 
(VDF) . Substituting /i = cos 6 yields 



V T - V r = Vfl, V t = Vg + Vj, = V (1 - fi ) (31) 

These coordinates are appropriate for a series expansion 
of the VDF in Legendre polynomials ( Larson|1970 1: 



f(V, tl )=g(V) + J2a l (V)Pi( f ,) 

oo 



(32) 



where A (V) = g(V) + a (V) and A,(V) = ai(V) for I > 1. 
In this expansion 



=-47rG 2 m./lnA 



dh 

dvi 



1 3 

E 



2 ^— ' dvidvj 

1,3 = 1 



d 2 9 
dvidvj 



(29) 



To compute this expression we need to know the phase 
space distribution function f(f, v, t). We approximate the 
distribution function by a series expansion which accounts 



for the spherical symmetry of the system (see e.g. Rosen- 
|bluth et al.][T9^7||Larson|[T970l ). The expansion coefficients 
are expressed in terms of the velocity moments needed to 
compute the collisional terms of the moment equations. We 
then compute the phase space distribution function / to 
calculate the Rosenbluth potentials h and g, the right-hand 
side of the FP equation and thus the collisional terms of the 
Boltzmann equation. 



4-2.1 Construction of the distribution function 

The phase space distribution function in spherical symmetry 



only depends on 



2 i 2 



and t, i.e. / = f(r,v r ,v e + 



u 0i 



t), which implies that the system is axially symmetric 
in the velocity space with axes v r , v$ and v^. The velocity 
components can be written in spherical coordinates 



1 V 2 
9(V) = P j= exp( -) 



Thus g(V)V 2 is the Maxwell-Boltzmann (MB) VDF. 
Pi(fJ,) are the Legendre polynomials, and the functions a,(V) 
are defined by 



Oi(V)= S (V) ^ 

2=0 

where l ma x denotes the highest order of the Legendre poly- 
nomials Pi(/z) in the expansion of the VDF. Due to axial 
symmetry in velocity space the VDF can only depend on 
powers of Vt and v r . Using equations ( 31 1 and fully expand- 



ing the VDF we find the following constraints for the coef- 
ficients c nm : 

• n ^ to 

• n and m are either both even or both odd 



We obtain for model b a VDF which extends to order 
I = 5 in the Legendre Polynomials Piifi) which reads 
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We obtain for the VDF of model a the coefficients ct 



f(V, n) = g(V) + g(V)(c 00 + c 02 V 2 + c 04 V 4 )P Q (p) 
+ g(V)( Cll V + c 13 V 3 + c 15 V 5 )Pi( M ) 

+ g(V)(c 2 2V 2 +C 2 4V 4 )P 2 (fl) 

+ g(V)(c 33 V 3 + c 35 V 5 )P 3 (fi) + g(V)c i4 V 4 P 4 (V) 
+ g(V)V 5 c 55 P 5 (fi) 

(33) 



The VDF for model a only extends to order I = 4 and can 
be obtained from equation ( 33 1 by setting all coefficients dj 
with j > 4 to zero. 

We can now calculate the coefficients Cij using the defi- 



nition of the centralized moments from equation 1 171. How- 



ever, we first have to transform equation (171 to the new 
coordinate system (V,p). The volume element d 3 v in these 
coordinates is written as: 



coo 

C02 
C04 



C24 



C33 



27 7(p r +2pt) K r + 2/trt + Kt 



C22 = 



4p<r 2 



8pa 4 



7 p r + 2p t K r + 2/trt + Kt 



4(7 2 



+ 



per- 



12po- 6 



1 p r + 2p t K r + 2/t r t + K t 



Cll = - 



C13 



8<7 4 12pcr 6 ' 120pcr 8 
F r + F t 

2pcr 4 
F r + Ft 

10pa 6 

3(Pr — Pt) 2k t + Krt — K t 



(38) 



2/9CT 4 



12pa 6 



Pr - Pt 2«r + Krt — Kt 



6pa 6 

- \Ft 
15 per 6 



84pa 8 



C44 = 



35pa 8 



d 3 u = y 2 dVd(cose)d<^ = V 2 dV d/x d(?i> where p = cos8 

(34) 

Thus we obtain for the centralized moments: 

{n, m) — J d 3 vf (v r — it) n Ui" 

V 2 dV dp dtj> f (v r — uY'v™ 
2tv j 1/ 2 dl/dp/V n p n (V 2 (l-p 2 )) m/2 
2tt [ dl/dp/V 2+n+m p n (l-^ 2 ) m/2 (35) 



Since the coefficients dj with i = 0, 1 only depend on 
sums of moments, we can find a definition for total moments 
(see section 



5.11 



The role of anisotropy comes into the open when going 
to higher-order coefficients, like C2j oc (p r — pt) oc a, where 
a = 1 — pt/pr is the anisotropy parameter. We envisage a 
system as isotropic in a "weak" sense if a = everywhere. 
Strong isotropy holds if the distribution has the strict de- 
pendence / = f(r, (v r - v r ) 2 + (vj + v 2 ) 2 ,t) = f(r, V, t) on 
the modulus of the velocity, which results in F r — Ft = for 
the radial and tangential energy flux, i.e. spherical symme- 
try in velocity space. We now compute the moments of fifth 
order with the VDF of model a. In this case the expansion 
in Legendre polynomials Pi (p) expands up to I — 4 and the 
fifth order moments are 



We can obtain a linear system of equations to be solved 
for the coefficients dj by computing the different moments 
via this equation with the expansion of the VDF from equa- 



tion (33 1. It must be noted that: 



• The first centralized moment vanishes, since (v r — v r ) 
v r — v r = 0, i.e. 



G r = 10cr 2 Fr, Grt = 2q 2 F r + 3cr 2 F t , Gt = Sa 2 Ft 

(39) 

As we saw before, the system of differential moment equa- 



tions ( 20 1 and ( 21 1 combined with the mass relation ( 25 1 was 



not complete. We can now close it by including the three re- 



lations in equation (39 1. In these equations the fifth-order 



moments G are expressed through lower-order moments. We 



now have a set of equations ((20 1, (21 1, (251 and (39l) that 



is numerically solvable. This set describes our model a. 
Combining the three relations we have 



(1,0) =0 



(36) 



15 

O 



(40) 



• Since there are two tangential directions we add a factor 
2 in the definition below 



2p« = (0,2> 



(37) 



We will see that the left-hand side of equation (140 1) appears 



in the coefficient C55 when the coefficients dj of the VDF of 
model b are computed. It then becomes clear that equation 
( |40[ ) is a result of setting C55 = since C44 is the highest 
coefficient of the VDF of model a. 
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For the VDF of model b we find 

27 _ 

= 8 

7 



coo 



7(p r + 2p t ) (n r + 2n rt + Kt) 



C02 = — 



4a 2 



+ 



4pa 2 
(p r + 2p t ) 
pa 4 



8pa 4 

(fi r + 2n r t + Kt) 

12 pa 6 



CQ4 = 



1 

80^ 



Cll 



Cl3 = 



Cl5 



C22 



C24 = 



9(F r 



(p r + 2pt) , («r + 2K rt + K t ) 

12pcr 6 



4pcr 4 
4(F r + F f ) 

5p<r 6 
_ F r + Fj 
20pa 8 
3(p r -p t ) 
2pa 4 

(Pr - Pt) 
6pa 6 



+ 



+ 120p<r 8 

G r + 2G r t + Gt 



+ 



8pa 6 

G r + 2G r t + Gt 

20pa s 

G r + 2G r t + Gf 

280pcr 10 

2k v Krt 



Kt) 



(41) 



+ 



12p<7 6 
84pcr 8 



fit) 



C33 



ll(F r 



Ft 



G r 



hG r 



30pa 6 



C35 = 



C44 



30pcr 8 



+ 



2 

30p<r 8 

- ^G r t §G 

270pcr 10 



■ Krt + 3 



C55 



35pa 8 

- 5G r t 

945pcr 10 



We have the same dependencies of the coefficients Cij 
on the sum of moments and relations that determine the 



degree of anisotropy, such as p r — pt, F r — f^t. As we pre- 



dicted before, we can obtain relation ( 40 1 by setting C55 = 0. 
Similarly, the relation |kv — K r t + |ftt = obtained by cal- 
culating the fourth order moments with the VDF used in 



Spurzem & Takahashi ( 1995 1 can be found in the coefficient 
C44 of the VDF for I = 4 and I — 5 again. 

Computing the moments of order n + m = 6 leads to 
the four equations: 



H r = 15pu 6 — 45u 4 p r + 15a 2 n r 
H r t = dpcr 6 — 12a 4 p r — 6o- 4 p t + 2o 2 n r + 6<T 2 K r 
H tr = 8ptJ 6 — 8<r 4 p r — I6a 4 pt + 8a 2 Krt + <J 2 nt 
H t = 48pcr 6 - 144cr 4 pt + 18a 2 m 



(42) 



These equations are the closure relations for model b. 
The complete set of equations of model b consists therefore 
of equations pOl, plL p2l, pB) and Ipl. 



By means of equation (42 \ we also find the relation 



— H r — 4Hrt + 3Htr ~ T.^t = 0, 

15 6 



(43) 



where the left hand side of this equation appears in the 
coefficient Ce6 if we take the Legendre expansion of the VDF 
up to / = 6. 



5 WEAK ISOTROPY, TOTAL MOMENTS 
AND ROSENBLUTH POTENTIALS 

In this section we identify different degrees of isotropy. They 
are specified by anisotropy parameters that can be found in 



with the conducting gas sphere model of |Giersz fc Spurzem| 
( |1994[ ); Spurzem & Aarseth ( 1996 1. In these two studies the 



authors use a VDF of second order I = 2 in Legendre poly- 
nomials Pi(p) in order to compute the collisional terms of 
their model equations. 

f(V,p) = g{V)P {pt) + ^^g(V)P 2 (p) (44) 
2pcH 

As we explained before, the definition of weak isotropy is 



F, 



= Pt 
3 



F t , 



(45) 



This concept of isotropy includes second and third order 
moments. In the case of weak isotropy, the VDF becomes 
the MB distribution g(V) since -Po(p) = 1. To generalize 
the definition of weak isotropy we retake the MB VDF. This 
VDF describes thermal equilibrium and is defined as 



f(Y,li) = g(V) = p 



/2na 3 

We then compute the two moments of second order, 

2 

Pr = PO- 

2p t = 2pa 2 



(46) 



(47) 



The factor 2 in front of pt accounts for two tangential direc- 
tions. We recover the known isotropy condition by dividing 
the second equation by two and then subtracting the two 
resulting equations, 



Pr - Pt = 



(48) 



For a spherical symmetric stellar system this relation 
describes the highest degree of isotropy. We define the 
lowest-order anisotropy parameter as 



a p = p r — pt ■ 



(49) 



Thus, computing second order moments with a zeroth 
order VDF produces the two relations in equation (1471) . This 



can be used to derive the isotropy condition in equation ( 48 1 



which appears as an anisotropy parameter a p in the second 
order VDF of equation ( [44] ) . 

We can now recover an expression for the velocity dis- 
persion a by simply adding the two equations in (47 1 and 
solving for per 2 : 



2 Pr + 2p t 



(50) 



The random kinetic energy e is defined as e = (p r + 
2pt)/2; then, applying the isotropic condition p = p r = pt, 
we find that e = | p' = | pa 2 . This is the equipartition 
theorem for / = 3 degrees of freedom, which states that 
in thermal equilibrium at a temperature T every degree of 
freedom contains the same amount of average energy e, = 

|fc B r=|po- 2 . 

In order to find isotropy relations for higher-order mo- 
ments, we use a second- and fourth-order VDF. The fourth- 
order VDF was computed in the previous section. The 
second-order VDF is 
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/(Km) =g(V) + g(V)(c 00 + c 02 v 2 )P (p) 

+ g(V)cuVPiOi) + g(V)c 22 V 2 P 2 (p) 



(51) 



We determine the coefficients cy and then compute the 
fourth order moments re: 



6a 2 p r 



re r = — 3pa 4 
K rt — —2pa 4 + 2a 2 p r + 2a 2 pt 



(52) 



Kt 



OCT 4 + 16a 2 pt 



We can assume that these relations constrain our de- 
gree of anisotropy, since the information contained in higher- 
order moments can be expressed by lower-order moments. 
Similarly, we use now these relations to compute isotropy 
conditions that reappear as linear combinations of the re's 
in the coefficients c 22 , c 2 4 and C44 of the fourth-order VDF. 
These linear combinations should vanish in case of isotropy 
and thus can be identified as the anisotropy parameters of 
fourth order. 

For the linear combination of re's in the coefficient C44 



we directly find by inserting equation ( 52 1 



3re rt + - re f — 



(53) 



For the linear combination of the re's in C22 and C24 we obtain 
in the same way 

2re r + Krt — Kt = <4> p r =Pt (54) 



If the conditions for weak isotropy defined in equation ( 45 1 
hold as well, the fourth-order VDF only depends on the Leg- 
endre polynomials Po(p) and Pi(p) and sums of moments 
p, F and re. 

We thus conclude that linear combinations of moments 
in the coefficients dj for 2 ^ i describe anisotropy param- 
eters. The order of an anisotropy parameter is equal to the 
order of moments it consists of. The degree of isotropy is 
hence determined by the lowest-order anisotropy parame- 
ters that vanish. We therefore introduce a new definition 
for weak isotropy by requiring that all anisotropy parame- 
ters vanish, which corresponds to demanding that the VDF 
depends only on the Legendre polynomials Po{p) and Pi(p). 

5.1 Total moments 

We are now in position to define the total centralized mo- 
ments of the VDF, since this has been totally determined. 

(O = (((«r - V r ) 2 +vl+V 2 ) 2 } 

= {{p 2 V 2 + V 2 {l-p 2 )f) = {V n ) (55) 



2tt 



V n+2 f(jM,V) dp dV, 



since ve 



0. With the help of ( 55 I we calculate the 



even moments, which we define as p, re and H, 
p = {V 2 } =p r + 2p t 
re = (V 4 ) = re r + Krt + Kt 

H = (V 6 ) = 105(pa 6 - <7 4 (p r + 2p t )) + 21a 2 (re r + re rf + re t ) 
= H r + 3{H r t + H tr ) + H t 

(56) 



In the last line we have employed the relations given by ( 42 1 
With p and re and VDF up to order I = 4 we find 



3 re 
10 P ^ 



<"*>-'£(!'■ 



(57) 



(V 8 ) = a 



It 
3 p 



+ 3re 



which is independent of the uneven moments F r , Ft, G r , G r t 
and Gt- We thus define the uneven total moments 



F=-(F r + F t ) 

G = G r + 2G r t + Gt 



(58) 



The factor 1/2 in the definition of F is chosen in order to ob- 
tain consistency with the physical interpretation of F. This 
becomes clear when we add the two differential equations 
for the radial and tangential pressure p r and pt in equation 
(20 1, where we find that F — i(F r + Ft) corresponds to the 
radial flux of random kinetic energy. 

With these definitions our the coefficients of the VDF 
Cij for i — 0,1 now only depend on total moments. 

5.2 Rosenbluth Potentials 

After having calculated the expansion coefficients for the 



VDF f(V,p) in section 4.2.1 we now can calculate the 



Rosenbluth potentials, given by 

2w 1 V 

h(V,i*) = (m + mf) J!! ffr|x VfdVfdpfd^ 

0-10 
2tt 1 V 

g(y,n) = mj I J' J' f(V f ,ti f )W-v f \V?dV f d» f d<f> 
-1 

(59) 

So as to integrate for h(V, p) we can make a multi-pole ex- 
pansion, i.e. 

00 1 



1 

\v-Hf\ 



i = m = -i V > 



where 



u< = min(v, v') 
vy — max(v, v'), 
and the spherical harmonics are defined in the usual way. 



(60) 
(61) 



Yi. 



21 + 1 (1 m)\ lm0 
4tt (l + m)\ 1 V V " (62) 



Yi t . M (e,<t,) = {-i) m Y^, A (e,<j>), 

with the associated Legendre polynomials P; m (cos#). 

We use p = cos# and insert equation ( |60[ ) into the 
Rosenbluth potential h(V,p) of equation (59 1. After inte- 
grating over cj>, the associated Legendre polynomials are re- 
duced to Legendre polynomials Pi(p) and we hence can ap- 
ply the orthogonality relation 
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^ Pi(f*f)Pk(tJ.f)dfj,f = S kt 



(63) 



To compare our results with the lower-order estimation 
of |Spurzem fc T akahashi ( 1995]) , we adopt their notation for 
the integrals over V, 



In = 
Kn = 



[ V V?g(V f )dV f 
Jo 



V?g(V f )dV f 



(64) 



With a VDF of order I = 5 we obtain the Rosenbluth 
potential h(V, fi): 

h(V,n) 



47r(m + rrif) 
h 



+ Kr) (1 + coo) + (y + A- 3 )c 02 + (y+Ks) 



('(14 



5V 
/ ho 



ho 



^3 



ftM 



V9V 5 ' 9 



ft(/i)+ 
c 5 sft(^), 



(65) 

If we set all coefficients dj — with the exception of C22, 



we recover the Rosenbluth potential h(V, /i) from Giersz & 



Spurzem ( 1994 1. This confirms the correctness of our result. 



Moreover, we obtain the Rosenbluth potential h(V, /x) for 
order I — 4 by setting the coefficients C15 = C55 = 0. 
To calculate g( V, (x) we write 



Vf\ 



(\v-v f \) 2 _ ( V 2 + Vf — 2VVf cos x) 



(66) 



\V-Vf\ \V-Vf\ 

where \ ls the angle between the vectors v and vf. This can 
be rewritten in terms of the angles 0, (j>, 6j, <f)f with the 
general formula for Legendre polynomials 

p i (cos X )= £ (j - M)j P^\ cos e)P^(cose f )e- im ^-^ 
m=-i V t + I m !i- 

(67) 



Setting I — 1 we can substitute cos \ m equation ( 66 1 with 



Pi(cosx) from equation (67 land insert the result into the 
Rosenbluth potential g(V, fi) of equation ( |59| |, This leads to 
an expression for g(V, fi) depending on products of Legen- 
dre and associated Legendre polynomials. After carrying out 
the integration over <f> we use relations between the Legendre 
and associated Legendre polynomials that reduce the prod- 
ucts and enable us to apply the orthogonality relation (j63h 



We can now write the result in the notation of |Spurzem fc| 
Takahashi ( 1995 1 to verify that their lower order potential 



g(V, fj.) is contained in our result by using a VDF up to order 
I = 5; 



(Vh + T^h + ^j-Ki + A 3 ) (1 + coo) + 

/IV 2 \ 

[Vh + -^yh + —A 3 + A5JC02+ 

IV 2 \ ' 

VI 6 + + -yA 5 + K 7 )c 04 



PoW 

1, 



VKa)cu+ 



3 6 15V 2 

3 15V 2 

h , Js_ 
15V 35V 3 



+ ^V 3 A 3 



\ VK *)** 



iy 3 if 6 -^if ? )rr, 



A(m) 



'•22 



A ho —V i K 
15V 35V 3 35 3 



P2M 



h 
35V 2 



ho V* 
63V 4 63 1 



V 3 
35 



A3 )C33 + 



+ 



+ 



ho I12 V „ 
35V 2 63V 4 63 3 



J10 /12 V 
63V 3 99V 5 99 1 



A 5 )c 35 



A2 ii4 



V 



_ ^ 
35 

6 y4 
63 

r 5 



ftM 



A3 C44 



ft (a*) 



99V 4 143V 6 143 



Ai 



99 



A 3 



C55 ft (m) 



(68) 



We again can set all coefficients dj = with the exception 
of C22- We find that this leads to the second-order result of 
|Spurzem fc Takahashi| |1995), which corroborates the cor- 
rectness of our result for g(V, fj,). The fourth-order Rosen- 
bluth potential g( V, fi) can be recovered by setting the coef- 
ficients C15 = c 55 = 0. 

Eventually we carry out the integration over V for both 
Rosenbluth potentials h( V, [i) and g(V, n) which is needed 
for the further computation of the right-hand sides of the 
moment equations. 



6 COLLISION TERMS 

With the coordinates in velocity space V and \i the Fokker- 



1957) 



Planck equation (|29[) transforms to (from [Rosenbluth et al. 



14 Schneider, Amaro-Seoane & Spurzem 



i { 6f(y,n) 
r\ St 



V 2 dV\ n ,W dV 



(69) 



1 d ( ( 2 \ dh(V,p) \ 

2d 2 g(V,^y 



+ 2^W^( m ' l)V dV 2 

+ — — (f(vJ {1 ~^ 2d29{v '^ 
2V 2 9 M 2 V V v 2 dp 2 

(l-p 2 )dg(V,p) (l-p 2 )dg(V,p) 

t t r r M 



+ 



V 

1 d 2 

V 2 dVdp 

1 d 

2V 2 dV 
H dg(V, p) 
V dp 

i d 



dV 

f(V,p)(l-p 



v 2 

2 . fd 2 g(V,p) ldg{V,p) 



/(Km) 



avs/i v dp 
(i- M 2 )a 2 5 (K, M ) d 5 (v» 



v 



dV 



+ 2\/ 2 9 M ^ j V ^ 2 5 M 2 

2 » dg{V,n) , 2 (l-M 2 )9 2 ff(KM) 2 9<;(V» 



r 



V 2 



9/i 



where T = 47rG 2 m/lnA and In A is the Coulomb loga- 
rithm. To obtain the collision terms of the moment equations 



(201, (21 1 and (22 1, we multiply the FP equation with pow- 



ers of velocity components and integrate over velocity space 
' 8{n, m) ' 



St 



= d 



V I — j (V r - V r ) V t 

V / cnc 



= 2tv [ V 2 dVdp f 6 ^^ ) V n p"(V 2 (l - p 2 )) m/2 

•* \ / cnc 

= 2tt f dV dp ( Sfi J^ ) ) V 2+n+m p n (l - p 2 ) m/2 

" V / cnc 



(70) 



For a single-mass model rrif = m and some collisional terms 
must vanish. These are the particle density p, due to parti- 
cle/mass conservation, the collision term of the bulk velocity 
u (or pu since internal collisions do not disturb the motion 
of the barycenter) and the collisional term for the energy 
density defined as e = (jp r + 2p t )/2, due to energy conserva- 
tion; 



= 0, 



= 0, 



5p r 
St 



+ 2 



5p t 
St 



(71) 



= 0, 



as expected, which proves that our calculations are right. 
We define the anisotropy parameters, which appear in the 



a P = Pr — Pt 
a F = 2F r - 3F t 

a K l = 2K r + K r i — Kf 

a K 2 = 8/t r — 24K rt + 3Kt 
aoi = 2G r — G r t — 3Gt 
a G2 = 8G r - 40G rt + 15G t 
and use the total moments 

F=±(F r + F t ) 

K = K r + 2/Crf + K t 
G = G r + %G r t + Gt 



(72) 



(73) 



We then give the collisional terms for the two models a and 
b in appendix [X] 



7 THE VELOCITY DISTRIBUTION 
FUNCTION 

In this section we investigate the influence of moments and 
anisotropy parameters on the VDF. For that, we use a VDF 
with moments up to fifth order. We express the VDF with 
the total moments F,k and G and anisotropy parameters 
a p , a F , a re i,2 and aai,2, 



f(V,p)=g(V) 



, T ,,# 15 5V 2 



8a 4 



+ 



8pa 4 l2pa G 

nr J WF 
+ 9(y \-2p^ 

VG V 3 G 

+ 



V 2 K V 4 K 
+ 



+ 



+ 



120pcr 8 
8V 3 F 

5pa e ~ lOpcr 8 



V 5 F 



V 5 G 



8pa 6 20pa s 280pa 1 

■ n „f 3V 2 a p V 4 a p V 2 a Kl V 4 a Kl \ 
+ 9( ' V) {^ 2 --l^i--12p^ + 8^) P2M 

UV 3 a F V 5 a F V 3 a G i , V 5 a G i \ D , 



60p<r 6 



+ 



30p<r 8 540pa 
V 5 a G2 



10 



(74) 

In order to obtain the MB VDF in the case of thermal 
equilibrium, g(V)V 2 , we have to multiply f(V,p) with V 2 . 
In the figures, the V-axis denotes the modulus of the veloc- 
ity and the /x-axis the direction of the velocity vector. When 
p = 0, the radial velocity component is ty = pV = and 
the tangential velocity component is vt = \/V 2 (l — p 2 ) = V 
and vice- versa for p, = 1. The z-axis indicates the phase 
space probability density f(V, p). If not stated otherwise, we 
choose a = 10kms _1 pc" 3 and k = 150 000km 4 s" 4 pc" 3 . We 
normalize f(V, p) by setting the particle density p = lpc -3 , 
and then J /(V, p)V 2 dVdp = 1. We set to zero the values of 
the moments F and G and the anisotropy parameters a v , a F , 
a K i,2 and cigi,2. To emphasize the effects of moments and 
anisotropy parameters, we choose very high and low values 
for these quantities in some plots. This results in negative 



Higher order moment models of dense stellar systems 15 



values of the distribution function which is unphysical but 
reflects the polynomial ansatz of the truncated series expan- 
sion of the VDF. We explore the parameter space to analyze 
their influence on the VDF. 

In order not to clutter the figures with the values for 



the set of parameters listed in equations ( 721) and ( 73 I the 
values for the plots are collected in table [T The parameters 
that change the VDF with respect to the MB distribution 
are denoted in the each plot. 

Figure [2] displays two plots of the VDF f(V,p)V 2 . In 
the left plot k — and it is clear due to the shape that this 
is not the MB VDF for thermal equilibrium. To choose the 



right value for re compute it by means of equation ( 55 1 using 
the MB distribution g(V)V 2 this yields^ 

re = 15/9cr 4 (75) 

For given a and p this is the value we have to choose for 
re which is in our case re = I50 000km 4 s- 4 pc" 3 . Then we 
obtain the MB distribution as can be seen in the right plot 
of figure [2] where this value was used. 

The anisotropy a p describes the difference between the 
second order moments p r and pt which represent the radial 
and tangential pressure (or equivalently energy density) re- 
spectively. In thermal equilibrium we have p r = p t and thus 
a p = 0. The second order moments determine the width of 
the VDF given by the dispersion a. When the anisotropy 
a p < the tangential pressure exceeds the radial pressure. 
As a consequence, we observe in the left plot of figure [3] that 
for fj, — t 1 the number of particles decreases whereas for 
H — > the number of particles increases. Since p determines 
the fraction of the radial and tangential velocity component 
this physically means that we have more particles with cir- 
cular orbits when a p < 0. For a p > we have the opposite 
behavior. 

A very similar effect is caused by the fourth order 
anisotropy a Kl as can be seen in figure[7] It appears together 
with Pi(ji) and the same powers of V as a p . However, since 
a K \ and a p appear with a different sign they have opposite 
effects. Consequently, we can assume that the fourth order 
anisotropy a K \ is a correction of the second order anisotropy 

The same argument holds for the third order anisotropy 
hf (figure[5| and fifth order anisotropy aai (figure 10 1. Both 
appear as factors of the Legendre polynomial P-z(p) with the 
same powers of V, but different sign. Thus aci can be seen 
as a correction of the third order anisotropy ap. 



In equation ( 74 1 we observe that uneven moments ap- 



pear with uneven Legendre polynomials Pi(p). However, 
these Legendre polynomials vanish at p = and thus the 
VDF is independent on uneven moments at p = 0. In other 
words, since p = corresponds to stars that have a vanish- 
ing radial velocity component v r , the distribution of stars 
that move on circular orbits are not affected by third order 
moments. This effect can be seen in in figures [4| [5] [9} |10| and 

El 

When including third order moments the VDF depends 



This is also the reason why |Louis| ( |l990F defined k' = k/15 in 
his model. Then thermodynamical equilibrium or isotropy yields 
k 1 = pa 4 . Nevertheless, this collides with equations ) |56| , where 
the total moments were computed giving a more natural defini- 
tion. 



on the total moment F and the anisotropy ap, where F 
is related to the total energy flux and describes differ- 
ences between radial and tangential energy fluxes. Negative 
values of the total moment F result in an increase of the 
maximum of the VDF for p, — > 1 (figure |4j left plot). We 
thus find more stars with eccentric orbits for F < 0. When 
F > the maximum of the VDF shifts to higher veloci- 
ties V when fj, — ¥ 1 (figure |4j right plot). This means that 
positive F increases the radial velocity component v r , but 
leaves the tangential component constant. As a consequence 
it increases the eccentricity of orbits, but not the number of 
stars with eccentric orbits as it does for F < 0. Note that 
the two plots in figure [4] have different scaling in the z-axis 
with respect to each other in order to display the distinct 
effects for F < and F > 0. Whereas the maximum of the 
VDF changes for F < from w 0.01pc~ 3 (km/s)~ 3 at p = 
to « 0.015pc~ 3 (km/s)~ 3 at p = 1 in the left plot of figure 
[4] it roughly stays constant for all values of p in the right 
plot corresponding to F > 0. The effect of the third order 
anisotropy ap is displayed in figure [5] < increases the 
number of stars with velocities above the mean and direc- 
tions corresponding to p w 0.5 whereas ap > causes an 
inverse effect. 

The fourth order moments are related to the kurtosis 
of the velocity distribution which gives a measure of high 
velocity stars as compared to thermal equilibrium. In figure 
[6] the VDF is plotted for two different values of the total 
moment of fourth order re. In the left plot re is chosen to 
be smaller than 15pa 4 . The VDF increases at its mean at 
the expense of high and low velocities. This corresponds to 
a deficiency of stars with high or low velocities but more 
stars with velocities near the mean, compared to thermal 
equilibrium. In the right plot the value of re is chosen to 
be higher than 15pa 4 . The wing of the VDF towards high 
velocities becomes thicker whereas the maximum of the VDF 
is smaller when compared to thermal equilibrium. Here the 
number of high velocity stars increases at the expense of 
stars with intermediate velocities. 

Whereas the anisotropy a K \ should be viewed in com- 
bination with the second order anisotropy Clp £LS mentioned 
before, the anisotropy a R 2 gives a new characterization of 
anisotropy at fourth order, which is displayed in figure [8] 

The effect of the total moment G on the VDF is illus- 
trated in figure [9] To get a physical understanding of this 
quantity we compare it to the total moment of third order 
F. F is an uneven moment which was considered to denote 
the radial flux of random kinetic energy. The random kinetic 
energy density e was given by the second order moments as 
e = (p r + 2p t )/2. Thus the third order moment F is the 
corresponding flux quantity for the second order moment p. 
Equivalently we can relate the fifth order moment G and the 
fourth order moment re. Since re is related to the number of 
high velocity stars G can be considered as a measure for the 



flux of these stars. Again aci (figure 10 1 should be viewed 
in the context of the third order anisotropy whereas aa2 



(figure 10 1 determines a new type of anisotropy at fifth or- 
der. 

Thus, every moment and anisotropy parameter has its 
own affect on the VDF. They act on different velocity inter- 
vals and redistribute stars from distinct orbitals. If we only 
include moments up to third order into our model, as it has 
been done in previous studies, our VDF is strongly limited. 
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Table 1. The table gives an overview over the different values of the total moments and anisotropy parameters in each plot. The values 
for the density p = 1 pc -3 and the velocity dispersion a = 10 kms" 1 pc — 3 are constant over all plots and therefore do not appear in the 
tabic. 




Figure 2. left: VDF where a = lOkms - 1 pc — 3 , the remaining moments and anisotropy parameters are set zero, right: VDF where 
a = lOkms - 1 pc -3 , k = 15p<T 4 , the remaining moments and anisotropy parameters are set zero. The right plot shows the MB distribution 
in thermal equilibrium 



We are then not able to describe areas in velocity space as 
is possible with moments of order > 3. More precisely, we 
obtain a much more detailed description of the distribution 
of stars in velocity space for stars with high velocities and 
stars which have neither radial nor tangential orbits, i.e. 
< ft < 1. 



8 DISCUSSION 

In this work we develop two statistical moment models for 
dense stellar dynamical systems. They are closed either at 
fifth- or at sixth-order depending on the required accuracy. 
They describe in a self-consistent way (including Fokker- 
Planck relaxation terms) local deviations of the velocity dis- 
tribution function from the MB distribution. The descrip- 
tion of the velocity distribution function includes third- and 
fourth-order moments. Third-order moments represent en- 
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Figure 3. VDF with two different values for the anisotropy a p . left: Shows the effect of negative anisotropy a p on a MB distribution. 
right: Shows the effect of positive anisotropy a p on a MB distribution. 




Figure 4. VDF with two different values for the third order total moment F corresponding to energy flux, left: Shows the effect of 
negative F on a MB distribution, right: The effect of positive F on a MB distribution 



ergy fluxes equivalent to asymmetries of the velocity dis- 
tribution around its center. Fourth order moments denote 
deviations from the MB distribution at high velocities. This 
cannot be described by a velocity distribution that is fully 
determined by its first two moments such as a Gaussian, 
commonly used to fit observational data. Due to the larger 
number of moments of the velocity distribution, the two 
models we introduce have the potential to fit detailed star- 
star and integrated light observations of globular clusters or 
nuclear star clusters in detail. However, they still underly a 
number of approximations, such as assuming spherical sym- 
metry, equal stellar masses, the Fokker-Planck and the local 
approximation and they also require a system with a high 
number of stars or high star densities so as to justify the 
statistical treatment. As the model equations only account 
for two-body relaxation, other mechanisms that drive the 
evolution of a stellar system can be added as terms in the 
model equations later on (e.g. unequal stellar masses and 
stellar evolution). In this work we have focused on giving 
the first complete analytical derivation of the relevant high- 
order moment equations. 

One of our goals is also to improve previous models such 
as the AGM or the moment model of |Louis| ( |1990[ ). For that, 
we achieve a more accurate modeling by including a larger 



number of moments. As we explained in section[7] increasing 
the number of moments leads to both a more complex VDF 
and an increasing number of differential moment equations 
(section |4.1[ ). This argument applies to the AGM rather than 
Louis' model. 

We therefore can describe the state of the system in 
terms of its phase space distribution function more accu- 
rately. As explained previously, GCs are in dynamic equilib- 
rium but not in thermodynamic equilibrium. While a system 
in thermodynamic equilibrium can be represented by a VDF 
that is fully defined by its first two moments, the number 
of non-vanishing moments increases for a system which is 
not in thermodynamic equilibrium. In most cases, it is im- 
possible to exactly compute the VDF for a system that is 
not in thermodynamic equilibrium. In a stellar dynamical 
system such as GCs or NCs there are numerous mechanisms 
that force the system away from thermodynamic equilib- 
rium, raising the issue of when to truncate the moment se- 
ries. Mechanisms that affect the high end of the VDF, such 
as the evaporation of stars from a stellar system, close three- 
body encounters and mass-segregation, suggest that the in- 
clusion of fourth- and fifth- order moments are important. 
This is also fortified by observations such as in the findings 
of high velocity stars in the core of Milky Way GCs. The 
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Figure 5. VDF with two different values for the third order anisotropy ap. left: Shows the effect of negative anisotropy ap on the VDF. 
right: The effect of positive anisotropy ap on the VDF. 




Figure 6. VDF with two different values for the total moments k. left: VDF with lower value for k with respect to the MB distribution. 
right: VDF with higher value for re with respect to the MB distribution. 



AGM, as a 3rd order model, does not accurately describe 
these mechanisms. 

The correct computation of the collisional terms is still 
a major difficulty and the local approximation is applied and 
an ansatz for the VDF is used to handle this problem. Even 
so, there are evident improvements over the previous mod- 
els that stem from the use of a larger number of moments 
(regarding the AGM) and a self-consistent method for the 
computation of the collisional terms (as compared to Louis 
1990). In contrast to the previous models, the collision term 



of a moment equation of order n does not only depend on the 
corresponding nth order anisotropy parameter but instead 
exhibits more dependencies on anisotropy of the parame- 
ters and moments of almost all orders as well. This leads to 
further coupling between the different moments. In a com- 
parative study between the AGM and FP and iV-body mod- 
els, [Spurzem~^TakahashT] |l995| concluded that in a multi- 
mass model a significant fraction of small-angle encounters, 
which transfer energy from the heavy to the light stars in the 
core, cause the light stars to move radially outwards on elon- 
gated orbits. As a result, the energy taken from the heavy 
particles is quickly redistributed over a much larger volume 
than assumed by the local approximation. Even though the 
local approximation is still applied in this model, the en- 



ergy transfer due to collisions should be improved due to a 
stronger coupling of the moments. This will provide a bet- 
ter estimate for the impact of the local approximation of the 
evolution of the system. 

The choice of the closure relation is very important, in 
particular at lower orders. In the AGM the system of equa- 
tions is closed with the heat flux equation, which relates the 
energy flux to the velocity dispersion. It is not clear how 
well the heat conduction closure of the AGM works. It obvi- 
ously allows the model to handle heat transfer and there are 
certainly parallels to gas-dynamics in GCs, but the descrip- 
tion of energy transfer via the gas-dynamical heat conduc- 
tion equation might nevertheless be a too simple description 
of this process. The heat conduction closure and the third 
order differential moment equations seems to be two com- 
pletely different descriptions of a similar process. Even so, in 
the comparative study by |Louis fc Spurzem| |T9"9"Tj ) between 
the AGM and the model of Louis ( 1990 1 reasonable agree- 



ment in pre-core-collapse could be achieved by proper choice 
of the free parameters of the AGM. However, these values of 
the free parameters of the AGM are not in agreement with 
the values resulting from the comparative study by |Giersz] 
|fc Spurzem| ( |1994[ ), where the AGM was fitted and com- 
pared to FP and TV-body models. This indicates that Louis' 
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a Kl =-7.5 * 1 4 km 4 s" 4 pc" 3 a Kl =1 .35 * 1 5 km 4 s" 4 pc" 3 




Figure 7. VDF with two different values for the fourth order anisotropy a K \. left: The effect of negative anisotropy a K ± on a MB 
distribution, right: The effect of a positive anisotropy a K i on a MB distribution 



a 1c2 =-4.5 * 1 5 km 4 s" 4 pc" 3 a K2 =6.0 * 1 5 km 4 s" 4 pc" 3 




Figure 8. VDF with two different values for the fourth order anisotropy a K 2- left: The effect of negative anisotropy a K 2 on a MB 
distribution, right: The effect of of anisotropy a K 2 on a MB distribution 



model does not agree very well with FP and TV-body mod- 
els. Furthermore, it has to be considered that the parameter 
A determining the heat conduction in gaseous models is just 
a scaling factor in isotropic gaseous models. In anisotropic 
gaseous models A prescribes the relative speed of the two 
relevant processes - the decay of anisotropy and the heat 
flow between warm and cold regions. With growing A heat 
flows faster, so there is less time for gravitational encoun- 
ters to destroy anisotropy ( Louis Sz Spurzem||1991 1. In our 
model (as in the model of Louis ( 1990)) this free parameter 
is absent. 



The closure equation in the model of Louis ( 1990 1 is 



an algebraic relation between the flux velocities of even mo- 
ments k, p and p. It is based on the assumption that the 
flux velocities of moments of order 2k increase with k. 

The closure relations we use are basically a mathemat- 
ical formulation of the fact that our model cannot describe 
an arbitrary degree of anisotropy, since its description of a 
stellar system is bounded by the highest moment that it in- 
cludes. It also reflects the limits of variability of the VDF 
and, thus, is a very natural choice. The only uncertainty of 
this closure relation is the error due to the polynomial ansatz 
for the VDF. The closure limiting the VDF is derived from 
the VDF itself. It does not stem from any other constraint 



that is independent of the form of the VDF arising from 
the boundary conditions like spherical symmetry and the 
absence of rotation. Hence, this ansatz should not be seen 
as an additional approximation, but rather as a consistency 
relation. 



The model equations consist of the set of equations ( 20 1 



and ( 21 1 for model a where the right-hand sides are given 
by equations JAlh to (IA6I) and the set of equations pOl, 



( 21 1 and ( 22 \ for model b with the right-hand sides given in 



equations I AT I to (A15). Furthermore, we need the Poisson 
equation | |25[ ) and the closure relations ( |39| ) or \42\ to com- 
plete the model equations. In order to exclude errors in the 
computation of the collisional terms, several measures were 
taken. The higher order Rosenbluth potentials were com- 
pared to the second-order Rosenbluth potentials of |Giersz"l 
& Spurzem (19941 and showed exact agreement. Further- 



more, the collisional terms for the density, bulk velocity and 
energy density vanish as expected according to mass and 
energy conservation and the fact that internal collisions do 
not disturb the motion of the barycenter. 

Eventually, several arguments have been given that pre- 
dict improvements of the fourth and fifth order models de- 
veloped in this work in comparison with its predecessors 
but a final estimate of the gained accuracy can only be ob- 
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Figure 9. VDF with two different values for the total moment G. left: Shows the effect of negative G on a MB distribution, right: Shows 
the effect of positive G on a MB distribution 




Figure 10. VDF with two different values for the fifth order anisotropy clqi. left column: Shows the effect of negative anisotropy clqi 
on the VDF. right column: Shows the effect of positive anisotropy clqi on the VDF. 



tained by means of numerical simulations and subsequent 
comparison with other models. The next step will be to im- 
plement and test the model in a numerical code such as the 
anisotropic gaseous modeQ For that, the left-hand sides of 
the differential moment equations ( |20[ ), ( |21[ | and ( |22[ ) have 
to be discretized as in the appendix of [Amaro-Seoane et al.| 



(20041, and the collisional terms which form the right-hand 



sides of the moment equations have to be reformulated. They 
should be simplified and reordered to allow an effective im- 
plementation into a numerical code. 
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Figure 11. VDF with two different vaiues for the fifth order anisotropy aci- left: Shows the effect of negative anisotropy ac2 on the 
VDF. right: The effect of positive anisotropy ag2 on the VDF. 
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